III. Soving the diffusion problem with implicit Euler integration with multigrid & finite differences
In [254]:
using PyPlot
For now, we want to solve the diffusion problem on a rectangular domain. The reaction part of the equation we want to solve in the end is currently set to zero. We therefore want to find a function $u$ defined as follows.
\begin{align} u: [a, b] \times [c, d] \times ℝ^+ &\to ℝ^+ \\ (x, y, t) &\to u(x, y, t) \end{align}We therefore want to solve the following PDE:
$$ u_t = D Δ u = D (u_xx + u_yy) \\ u(a,y,t) = u(b,y,t) = u(x,c,t)=u(x,d,t)=0 $$For this fourth solution, we still use a second-order finite differences method in space. In time, however, we use an explicit singly diagonal implicit Runge-Kutta method. The details of this integration step are mentioned in the corresponding section below.
For the spatial derivatives, we don’t solve for $u^{n+1}$ directly (e.g. with $LU$-decomposition). Instead, we use a multigrid approach to solve the system iteratively. As the spatial approximation with finite differences has a simple, regular form, we can easily create the finite differences matrix for other, coarser grids as well. The matrices have the exact same form as for the finest grid, except that they have fewer entries. However, we impose that the number of intervals is a power of two to ensure that they can easily be mapped to a coarser grid.
In [255]:
immutable Grid2D
x_start::Real
x_end::Real
y_start::Real
y_end::Real
Nx::Int
Ny::Int
dx::Real
dy::Real
levels::Int
function Grid2D(x_start, x_end, y_start, y_end, Nx, Ny, levels)
if Nx/2^(levels-1) != int(Nx/2^(levels-1))
error("ERROR: The number of intervals in x-direction is not divisible by 2^i")
end
if Ny/2^(levels-1) != int(Ny/2^(levels-1))
error("ERROR: The number of intervals in y-direction is not divisible by 2^i")
end
new(x_start, x_end, y_start, y_end, Nx, Ny, (x_end-x_start)/Nx, (y_end-y_start)/Ny, levels)
end
end
In [256]:
immutable TimeSteps
t_start::Real
t_end::Real
Nt::Int
dt::Real
TimeSteps(t_start, t_end, Nt) = new(t_start, t_end, Nt, (t_end-t_start)/Nt)
end
As a very simple scenario, we use the following initial conditions:
$$ u(x,y,t = 0) = C_0 \cdot sin(π \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) $$This leads to the following analytical solution:
$$ u(x,y,t) = C_0 \cdot sin( \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) \cdot exp( -\left( \frac{1}{(b-a)^2} + \frac{1}{(d-c)^2} \right) \pi^2 D t) $$
In [257]:
function gen_exact(grid, C0, D)
x_size = grid.x_end - grid.x_start
y_size = grid.y_end - grid.y_start
exact(x,y,t) = C0 * sin((x-grid.x_start)/x_size * pi) * sin((y-grid.y_start)/y_size * pi) *
exp( - (1/x_size^2 + 1/y_size^2) * pi^2 * D * t)
end
Out[257]:
We define the finite differences matrices seperately for the $\partial_{xx}$ and $\partial_{yy}$ operators. The vector of unknowns $u$ is organized in col-major fashion (neighboring elements have the same y-value but not the same x-value).
Thus, the matrix for $\partial_{xx}$ is block-diagonal and organized as follows:
\begin{align*} FD_x = \begin{pmatrix} B_1 & 0 & \cdots & 0 \\ 0 & B_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & B_{N_y} \end{pmatrix} & \hspace{1cm} \text{with} & B_i = \begin{pmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & 1 & -2 \end{pmatrix}_{N_x \times N_x} \end{align*}The matrix for $\partial_{xx}$ has has non-zero entries in the diagonal and the $N_{xx}$-th upper and lower diagonals.
\begin{align*} FD_y = \begin{pmatrix} -2 & 0 & \cdots & 1 & \cdots & 0 \\ 0 & -2 & 0 & \ddots & \ddots & \vdots \\ \vdots & 0 & -2 & \ddots & \ddots & 1 \\ 1 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 1 & \cdots & 0 & -2 \end{pmatrix}_{N_x N_y \times N_x N_y} \end{align*}The sum of these matrices, combined with their step sizes, forms the finite difference matrix that is used for the integration.
$$ FD = \frac{1}{\Delta x^2} FD_x + \frac{1}{\Delta y^2} FD_y $$As we have to create several different matrices for the multigrid matrix, we define a function that returns those matrices depending on the number and size of intervals.
In [258]:
function fdmatrix(grid::Grid2D, matrix_function::Function = A -> A)
# number of unknowns in x and y direction
ux = (grid.Nx-1)
uy = (grid.Ny-1)
# finite difference matrix for ∂²u/∂x²
FDx_local = spdiagm((ones(ux-1),-2*ones(ux),ones(ux-1)),(-1,0,1))
FDx = blkdiag({FDx_local for i in 1:uy}...)
# finite difference matrix for ∂²u/∂y²
FDy = spdiagm((ones(ux*(uy-1)),-2*ones(ux*uy),ones(ux*(uy-1))),(-ux,0,ux))
# return combined matrix for Δu
matrix_function(1/grid.dx^2 * FDx + 1/grid.dy^2 * FDy)
end
Out[258]:
In [259]:
# helper function to get coarser levels of a grid
function level(grid::Grid2D, lvl)
if lvl > grid.levels
error("ERROR: Grid doesn't support level $(level).")
end
Grid2D(grid.x_start, grid.x_end, grid.y_start, grid.y_end, int(grid.Nx/2^(lvl-1)), int(grid.Ny/2^(lvl-1)), grid.levels-(lvl-1))
end
Out[259]:
In [260]:
function weighted_jacobi(A, x, b)
weight = 2/3
(speye(A) - weight*A) * x + weight*b
end
Out[260]:
In [261]:
function restrictionmatrix(grid::Grid2D)
ux_old = grid.Nx - 1
uy_old = grid.Ny - 1
ux_new = div(grid.Nx,2) - 1
uy_new = div(grid.Ny,2) - 1
R = spzeros(ux_new*uy_new, ux_old*uy_old)
for j = 1:uy_new
for i = 1:ux_new
i_new = ux_new*(j-1) + i
R[i_new, ux_old*(2*j-1) + 2*i ] = 4
R[i_new, ux_old*(2*j-1) + 2*i - 1] = 2
R[i_new, ux_old*(2*j-1) + 2*i + 1] = 2
R[i_new, ux_old*(2*j-1 - 1) + 2*i ] = 2
R[i_new, ux_old*(2*j-1 - 1) + 2*i - 1] = 1
R[i_new, ux_old*(2*j-1 - 1) + 2*i + 1] = 1
R[i_new, ux_old*(2*j-1 + 1) + 2*i ] = 2
R[i_new, ux_old*(2*j-1 + 1) + 2*i - 1] = 1
R[i_new, ux_old*(2*j-1 + 1) + 2*i + 1] = 1
end
end
R /= 16
end
Out[261]:
In [262]:
function gen_mgsolver(matrix_generator, restriction_generator, smoother, grid::Grid2D, ν1::Int, ν2::Int, tol=1e-12, max_it = 100,)
# create main matrices for all levels
A = Dict()
for i = 1:grid.levels
A[i] = matrix_generator(level(grid,i))
end
# create restriction matrices for all levels
R = Dict()
for i = 1:grid.levels-1
R[i] = restriction_generator(level(grid,i))
end
function mgsolver(b::Array, x0::Array = b, level::Int = 1)
if level == grid.levels
return A[level] \ b
end
x = x0
for it = 1:max_it
# smoothen problem on current level
for i = 1:ν1
x = smoother(A[level], x, b)
end
# restrict residual to next level and solve recursively
r = R[level] * (b - A[level] * x)
e = mgsolver(r, zeros(r), level+1)
# prolongate error to current level and add to solution
x += 4*R[level].' * e
# post-smoothening on current level
for i = 1:ν2
x = smoother(A[level], x, b)
end
if level > 1
return x
end
if norm(A[level] * x - b) < tol
#println("Converged with $(it) iterations")
return x
end
end
return x
end
return mgsolver
end
Out[262]:
We define an ESDIRK method with the following coefficients. As the implicit coefficient σ is the same for all stages, we always have the same matrix in the system $Ax=b$ that we solve – the only thing that changes is the right-hand side. To simplify the formulation, we use the following notation:
$$ FD = \text{finite differences matrix} \\ A = \Delta t \cdot D \cdot FD $$This way, $\Delta t \cdot f(u) = A \cdot u$ and we get the following formulation for the intermediate stages $k_i = \Delta t \cdot f(\xi_i)$:
$$ k_1 = A \cdot u^n \\ k_2 = A \cdot \left( I - \sigma A \right)^{-1} \left( u^n + a_{21} \cdot k_1 \right) \\ k_3 = A \cdot \left( I - \sigma A \right)^{-1} \left( u^n + a_{31} \cdot k_1 + a_{32} \cdot k_2 \right) \\ k_4 = A \cdot \left( I - \sigma A \right)^{-1} \left( u^n + a_{41} \cdot k_1 + a_{42} \cdot k_2 + a_{43} \cdot k_3 \right) $$Finally, we can compute the next time step of $u$:
$$ u^{n+1} = u^n + b_1 \cdot k_1 + b_2 \cdot k_2 + b_3 \cdot k_3 + b_4 \cdot k_4 $$
In [263]:
function gen_esdirk3(D::Real, grid::Grid2D, time::TimeSteps, matrix_generator::Function, solver::Function)
# generate matrix for RHS: f(u) = Au
A = time.dt * D * matrix_generator(grid)
# coefficients for ESDIRK3
c = [0 1767732205903/2027836641118 3/5 1] # not used
σ = 1767732205903/4055673282236
a1 = [0 0 0 0]
a2 = [σ σ 0 0]
a3 = [2746238789719/10658868560708 -640167445237/6845629431997 σ 0]
a4 = [1471266399579/7840856788654 -4482444167858/7529755066697 11266239266428/11593286722821 σ]
b = a4
#=
order_conditions = [
sum(b) - 1
b * c.' - 1/2
b * (c.^2).' / 2 - 1/6
b[1] * sum(a1 .* c) + b[2] * sum(a2 .* c) + b[3] * sum(a3 .* c) + b[4] * sum(a4 .* c) - 1/6
]
println(order_conditions)
=#
function esdirk3(u)
k1 = A * u
k2 = A * solver(u + σ * k1)
k3 = A * solver(u + a3[1] * k1 + a3[2] * k2)
k4 = A * solver(u + a4[1] * k1 + a4[2] * k2 + a4[3] * k3)
u += b[1] * k1 + b[2] * k2 + b[3] * k3 + b[4] * k4
end
return esdirk3
end
Out[263]:
In [264]:
function integrate(x::Array, integration_step::Function, time::TimeSteps)
for i in 1:time.Nt
x = integration_step(x)
end
x
end
Out[264]:
In [265]:
g = Grid2D(1,3,2,5,2^5,2^5,3)
Out[265]:
In [266]:
t = TimeSteps(0, 1, 200)
Out[266]:
In [267]:
C0 = 5
D = 0.1
exact = gen_exact(g, C0, D)
Out[267]:
In [268]:
C_exact = Float64[exact(x,y,t.t_end) for x=g.x_start:g.dx:g.x_end, y=g.y_start:g.dy:g.y_end]
pcolormesh(C_exact)
colorbar()
clim((0,5))
In [269]:
# use exact function to set up initial conditions
C = Float64[exact(x,y,t.t_start) for x=g.x_start:g.dx:g.x_end, y=g.y_start:g.dy:g.y_end];
In [270]:
σ = 1767732205903/4055673282236
solver = gen_mgsolver(g -> speye((g.Nx-1)*(g.Ny-1)) - t.dt * σ * D * fdmatrix(g),
restrictionmatrix, weighted_jacobi, g, 0, 1)
integration_step = gen_esdirk3(D, g, t, fdmatrix, solver)
Out[270]:
In [271]:
x = vcat(C[2:end-1,2:end-1]...)
x = integrate(x, integration_step, t)
C[2:end-1,2:end-1] = reshape(x,g.Nx-1,g.Ny-1)
pcolormesh(C)
colorbar()
clim((0,C0))
In [272]:
norm(C - C_exact)
Out[272]: